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Several new families of nonlinear three-dimensional travelling wave solutions to the 
Navier-Stokes equation, also known as exact coherent states, are computed for New¬ 
tonian plane Poiseuille flow. The symmetries and streak/vortex structures are reported 
and their possible connections to critical layer dynamics examined. While some of the 
solutions clearly display fluctuations that are localized around the critical layer (the sur¬ 
face on which the streamwise velocity matches the wave speed of the solution), for others 
this connection is not as clear. Dynamical trajectories along unstable directions of the 
solutions are computed. Over certain ranges of Reynolds number, two solution families 
are shown to lie on the basin boundary between laminar and turbulent flow. Direct com¬ 
parison of nonlinear travelling wave solutions to turbulent flow in the same channel is 
presented. The state-space dynamics of the turbulent flow are organized around one of 
the newly-identified travelling wave families, and in particular the lower branch solutions 
of this family are closely approached during transient excursions away from the dom¬ 
inant behaviour. These observations provide a firm dynamical-systems foundation for 
prior observations that minimal channel turbulence displays time intervals of “active” 

turbulence punctuated by brief periods of “hibernation” (see e.g. Xi, L. and Graham, 
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M. D., Phys. Rev. Lett., 104 , 218301 (2010)). The hibernating intervals are approaches to 
lower branch nonlinear travelling waves. Representing these solutions on a Prandtl-von 
Karman plot illustrates how their bulk flow properties are related to those of Newtonian 
turbulence as well as the universal asymptotic state called maximum drag reduction 
(MDR) found in viscoelastic turbulent flow. In particular, the lower and upper branch 
solutions of the family around which the minimal channel dynamics are organized appear 
to approach the MDR asymptote and the classical Newtonian result, respectively, both 
in terms of bulk velocity and mean velocity profile. 

1. Introduction 

The understanding of the nature of near-wall turbulence has been greatly advanced by 
recent applications of dynamical systems theory to turbulent flow (Kawahara et al. 2012). 
In particular, over the past two decades, the discovery of three-dimensional fully nonlinear 
travelling wave (TW) solutions to the Navier-Stokes equations has enabled a priori study 
of self-sustained near-wall coherent structures that resemble in many ways that transient 
structures observed in fully turbulent flows (Hof et al. 2004). These solutions, also denoted 
as exact coherent states (Waleffe 2001) (ECS), are steady states in a reference frame 
translating at a constant streamwise speed. They have been found numerically in all 
canonical wall-bounded geometries for turbulent flows (plane Couette and Poiseuille, 
pipe and boundary layer) (Nagata 1990; Clever & Busse 1997; Nagata 1997; Waleffe 
1998, 2001, 2003; Faisst & Eckhardt 2003; Wedin & Kerswell 2004; Gibson et al. 2008, 
2009; Schneider et al. 2010; Duguet et al. 2012; Blackburn et al. 2013). Most solutions 
that have been found to date are spatially extended, but recent studies show the existence 
of spatially-localized travelling solutions that closely resemble turbulent puffs in the pipe 
flow geometry (Avila et al. 2013; Chantry et al. 2014) or turbulent spots in the plane 
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Couette and Poiseuille geometries (Tillmark & Alfredsson 1992; Barkley & Tuckerman 

2005; Lemoult et al. 2013; Brand & Gibson 2014). Other exact coherent states or other 

types of invariant solutions to the governing equation have been also found numerically. 

Periodic or relative periodic orbits were computed for plane Couette (Kawahara & Kida 

2001; Viswanath 2007; Cvitanovic & Gibson 2010) and Poiseuille flow (Toh & Itano 

2003), pipe flow (Duguet et al. 2008) and two-dimensional Kolmogorov flow (Chandler 

& Kerswell 2013). Taken together, these solutions seem to form at least in part the 

dynamical skeleton underlying the chaotic dynamics of turbulent flow (Gibson et al. 2008; 

Kawahara et al. 2012). In the present work, we report and analyze several new families 

of such solutions in the plane Poiseuille geometry and further develop the understanding 

of the relationship between these solutions and the dynamics of turbulence. 

We focus here on plane Poiseuille (channel) flow of a Newtonian fluid with dynamic 
viscosity /r, density p and kinematic viscosity v = p,/pin a channel of half-height h. Char¬ 
acteristic inner scales are the friction velocity Ur = {t^I pfll"^ and the near-wall length 
scale or wall unit = vjur, where is the time- and area-averaged wall shear stress. 
As usual, quantities nondimensionalized by these scales are denoted with a superscript 
The friction Reynolds number is defined as Rcr = Urhjv = h/Si,. 

Exact coherent states primarily arise in pairs via a saddle-node bifurcation at a par¬ 
ticular Reynolds number. At the bifurcation, the pair of solutions emerges at finite am¬ 
plitude; we refer to each such pair as a solution family. Figure 1 illustrates a bifurcation 
diagram of solution amplitude versus Reynolds number for one such family (the “P4” 
family described below), using the maximum over y of the root mean square wall-normal 
velocity fluctuations normalized by the friction velocity as a measure of solution 

amplitude. An overbar indicates averaging over the streamwise and spanwise directions. 
The so-called lower branch (LB) solution of each pair denotes a low-drag state owing to 
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Figure 1. (Colour online) A bifurcation diagram for one solution family of travelling waves 
found in the present study (labelled P4 below), where the maximum in the root mean square 
wall-normal fluctuations is shown against Reynolds number. 

its lower maximum wall-normal velocity fluctuation compared to its corresponding up¬ 
per branch (UB) solution. (Additional solution branches can and do bifurcate off these 
primary states - the “P2” solutions described below are one such example.) In general, 
these solutions have a spatial structure in the form of low-speed streaks that are wavy 
in the streamwise direction, straddled by counter-rotating streamwise-aligned vortices: 
that is, they have the same basic qualitative structure as near-wall turbulence. 

The basic self-sustaining process underlying these structures has been qualitatively 
described by Waleffe (1997). More recently it has been observed for Couette flow that 
at least one lower branch solution family has a structure that consists of streaks, rolls 
and a weak streamwise-varying wave that develops a critical layer - i.e. its structure is 
localized around the surface where the streamwise velocity equals the wave speed of the 
ECS (Wang et al. 2007). In the classical linear stability theory of parallel shear flows, a 
critical layer is a planar surface around which normal mode perturbations localize (Drazin 
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& Reid 1981), while here the critical layer is a curved surface in 3D. Wang et al. (2007) 

presented a scaling analysis suggesting that the wavy fluctuations should be localized 

in a region of thickness 0{Re~^/^) and showed that this scaling was followed by their 

numerical solutions. In fact, they found that the flow structures at Re = 50000 and 

Re = 3000 were virtually identical modulo a Re^^^ rescaling of the direction normal to 

the critical layer surface. Hall and coworkers (Hall & Sherwin 2010; Blackburn et al. 2013) 

used a mixture of asymptotics and numerics to show, again for Couette flow, that the 

critical layer fluctuations couple back to the streamwise rolls to generate the nonlinear 

self-sustaining process that supports exact coherent states. In their formulation, this 

process is a version of wave-vortex interaction. They note that “remarkable” agreement 

is obtained between the high Re asymptotics and the numerical results down to Reynolds 

numbers of order 10^ (Hall & Sherwin 2010) - we emphasize this point because this is 

the Reynolds number range of the present results. Other recent and interesting work on 

critical layers and ECS is found in Viswanath (2009), Deguchi et al. (2013), Gibson & 

Brand (2014) and Deguchi & Hall (2014). Because of the clear importance of critical 

layer dynamics for at least some families of nonlinear travelling waves even at low Re, 

section 3.2 of the present work focuses on this topic. 

For plane Poiseuille flow, the first known families of travelling wave solutions were 
obtained by a homotopy continuation from ECS in plane Couette flow. Solutions were 
sought in a travelling reference frame using a Newton-Raphson method with the wave 
speed as an unknown (Waleffe 1998, 2001, 2003; Nagata & Deguchi 2013; Cibson & 
Brand 2014). Solutions have been also computed by a multiple shooting method (Itano 
& Toh 2001). It is worth noting that the lowest bifurcation point or onset Reynolds 
number for these solutions. Re ~ 977 or Rer ~ 44.3, is in good quantitative agreement 
with the Reynolds number for transition to turbulence observed in an experiment in this 
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geometry (Carlson et al. 1982). In fact, the solutions survive slightly below the critical 
Reynolds number for turbulence onset. Furthermore, the optimal spanwise wavelength 
105.5(5,^ at the onset of travelling wave solutions is remarkably matched well with that of 
experimentally and numerically observed near-wall streak spacing of 80 — 1205,^ (Smith 
& Metzler 1983). In the present work, several more families of travelling waves are found, 
which seem to have a closer connection to turbulent dynamics than the ones found earlier. 

An important issue regarding exact coherent states is their connection to the laminar- 
turbulent boundary - the boundary in state space between the basins of attraction of the 
laminar and turbulent states. Some of the lower branch ECS found in turbulent shear 
flows are embedded on this boundary (Itano & Toh 2001; Skufca et al. 2006; Wang et al. 
2007; Schneider et al. 2007; Kerswell & Tutty 2007; Duguet et al. 2008; Viswanath & 
Cvitanovic 2009). In particular, such solutions that have only one unstable direction are 
called edge states (Skufca et al. 2006). Because such solutions are somehow the weakest, 
most marginal form of self-sustaining turbulence, the structure of the basin boundary and 
dynamical trajectories that lie on it are likely to play an important role in understanding 
the dynamics of transition to turbulence or onset of turbulence in wall-bounded shear 
flows. Recently, an experimental observation has been reported for the existence of edge 
states in pipe flow (de Lozar et al. 2012). 

Returning to the dynamical-systems point of view, exact coherent states are periodic 
(or more complicated but still invariant) orbits in state space, while the time evolution 
of a turbulent flow is a dynamical trajectory wandering around them. An important 
question is how closely the turbulent trajectories approach these invariant states. For 
plane Couette flow, Gibson et al. (2008) visualized a clear illustration of this dynamical- 
systems viewpoint of turbulence by projecting the trajectory onto a set of orthonormal 
basis-states constructed with earlier ECS. Kerswell & Tutty (2007) also showed a clear 
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visual illustration for a pipe flow. In addition to the state space visualization, they pro¬ 
posed quantitative measurements of the distance between a given instant on the turbulent 
trajectory and ECS, and suggested that ECS are visited for approximately 10 % of the 
time in turbulent pipe flow. For a channel flow, the connections have yet to be fully made 
and will be investigated in the present study. 

One important motivation for gaining a better understanding of turbulence is the 
possibility of reducing drag. In this context, the lower branch ECS are attractive due to 
their low-drag flow features and a natural question is whether it might be possible to able 
to somehow steer turbulent trajectories toward these states. One very successful approach 
to turbulent drag reduction is to add small amounts of rheologically active additives such 
as flexible long-chain polymers into a liquid (White & Mungal 2008; Graham 2014). The 
most dramatic effect of the polymer additives on turbulence occurs in the near-wall 
region, weakening the turbulent eddies in this region. The key feature of these polymer 
solutions in drag reduction is an existence of the so-called maximum drag reduction 
(MDR) phenomenon, at which very high levels of drag reduction are achieved by polymer 
additives, first identified by Virk (Virk 1975). The most intriguing observation for MDR 
is its universal mean velocity profile, the experimentally observed upper limit on the 
amount of drag reduction that can be achieved with polymer additives, also known 
as the Virk asymptote (Virk 1975). This asymptotic limit is insensitive to changes in 
the polymer solution such as concentration, molecular weight or polymer type. Thus, 
for a given situation, the maximum amount of drag reduction achievable with polymer 
additives is invariant. 

Li et al. (Li et al. 2006; Li & Graham 2007) have investigated the effects of polymer 
additives on the channel flow ECS discovered by Waleffe (2001, 2003). For this solution 
family, as the level of viscoelasticity is increased, the Reynolds number for the solutions 
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to come into existence increases. The primary effect of viscoelasticity on these ECS is 
the weakening of the streamwise vortices. Other effects are also seen in changes in the 
budgets of turbulent kinetic energy, Reynolds stress and the mean shear stress. All these 
effects show, at least at low levels of drag reduction, that the basic mechanism of drag 
reduction by polymers can be clearly elucidated by examining the impact of polymers 
on travelling wave solutions. Nevertheless, these studies were limited to a single solution 
family and to relatively low Reynolds numbers and levels of viscoelasticity. 

Another set of recent studies, while not directly focused on ECS, sheds some light 
on the state-space dynamics of Newtonian and viscoelastic channel flow. Xi & Graham 
(2010a, 2012&) performed DNS studies of minimal channel flow at low Reynolds num¬ 
bers, finding that in Newtonian flow and for low to intermediate values of Weissenberg 
number, the flow cycles stochastically between “active” intervals, with strong streamwise 
vortices and three-dimensionality, velocity and a velocity profile near the von-Karman 
profile, and “hibernating intervals”, with very weak turbulence and a mean velocity pro¬ 
file approaching the Virk MDR profile. In viscoelastic flow, the polymers stretch during 
the active intervals, working against the streamwise vortices and shortening the duration 
of these intervals, while during the hibernating intervals the flow kinematics are very gen¬ 
tle and the polymers relax, only to begin stretching again at the beginning of the next 
active interval. Thus as the degree of viscoelasticity (Weissenberg number) increases, 
the overall dynamics look increasingly “Virk-like” as the active intervals contribute a 
decreasing time duration to the overall statistics. A simple theory is developed, based 
on exponential stretching of polymers during active intervals and the idea that these 
intervals cannot persist once the polymer stress reaches a threshold value. This theory 
predicts the Weissenberg-number dependence of the duration of the active intervals in 
good agreement with the simulations. At high Weissenberg number, the hibernation in- 
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tervals themselves are substantially altered and stabilized by viscoelasticity (Wang et al. 

2014), through mechanisms that are not yet understood. In addition, at low Reynolds 

number in the minimal channel flow geometry, hibernating turbulence is closely related 

to an edge state in which the Virk profile also arises, not just as a transient, but in the 

time-averaged velocity (Xi & Graham 2012a), even in Newtonian flow. A comprehensive 

overview of these studies is provided in Graham (2014). The connection between active 

and hibernating turbulence and upper and lower branch Newtonian EGS will be solidified 

in the present work. 

There are other indications as well that the Virk asymptote is not just universal for 
drag reduction by polymers but also arises in Newtonian turbulence. Dubief et al. (2011) 
observed in a simulation of Newtonian boundary layer flow that at a spatial position just 
upstream of where vortices and turbulence spots form, the mean velocity profile looks 
strikingly similar to the Virk MDR profile. Furthermore, the Virk MDR profile is also 
observed in a smoothed version of Newtonian plane Poiseuille flow in which spanwise 
length scales of the flow field below a specified size are suppressed (Kerswell et al. 2003). 
Finally, experimental observations of a Newtonian turbulent boundary layer flow sub¬ 
jected to spanwise wall oscillations display a mean velocity profile that, for y~^ < 30, 
closely resembles the Virk MDR profile (Bandyopadhyay 2006). 

In this paper, we present five new families of nonlinear travelling wave solutions in 
Newtonian plane Poiseuille flow, and examine their spatiotemporal structure and con¬ 
nections to the dynamics of turbulent flow in the same geometry. In particular, we find a 
family whose upper and lower branch solutions have mean velocity profiles that resemble 
Newtonian turbulent (von Karman) and MDR (Virk) profiles and we show the relation¬ 
ship between those solutions and trajectories of turbulent flows. The problem formulation 
and solution methodologies are presented in Section 2. Section 3.1 presents an overview 
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of mean flow properties and spatial structures of the solutions, while section 3.2 illus¬ 
trates the relation between the new travelling waves and critical layer dynamics. Sections 
3.3 and 3.4 describe the connections between the travelling waves, the laminar-turbulent 
boundary and turbulent dynamics. Section 4 presents conclusions. 

2. Formulation and solution approach 

We consider an incompressible Newtonian fluid in the plane Poiseuille geometry, driven 
with a constant mass flux Q. The characteristic length and velocity scales are the half¬ 
channel height h and the laminar centerline velocity Uc = (3/4)Q/h for the same mass 
flux, respectively. With these characteristic scales, the Navier-Stokes equations in nondi- 
mensional form are 


( 2 . 1 ) 


V • tt = 0 


du ^ 1 ^2 

——I- u ■ Vit = — Vp -I- --— V^tt. 
ot KCc 


( 2 . 2 ) 


Here, we define the laminar equivalent Reynolds number for the given mass flux as 
Rcc = Uchfv. Note that Ret = Uth/v = |i?ec, where Ub is the bulk velocity. The x,y 
and z coordinates are aligned with the streamwise, wall-normal, and spanwise directions, 
respectively. Periodic boundary conditions are imposed in the x and z directions with 
fundamental periods and and no-slip conditions are imposed at the walls y = ±1. 
The computational domain is thus [0,Lj;] x [—1,1] x [0,Lz] or simply [Lx,2,Lz\. The 
velocities are m, v, and w in the a;, y, and z directions, and the velocity at point (x, y, z) 
and time t is expressed as u = [it, v, ii'](a;, y, z, t). 

Computation of nonlinear travelling waves is performed using the open source code 
ChannelFlow written by Gibson (2012), with use of a Newton-Krylov-hookstep algorithm 
(Viswanath 2007). A numerical grid system is generated on x Ny x Nz (in x, y, and z) 
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meshes, where a Fourier-Chebyshev-Fourier spectral spatial discretization is applied to 

all variables. A typical resolution used is {Nx, Ny, N^) = (48,81,48). A travelling wave 

solution has the following form: 

u{x,y,z,t) = u{x - Cxt,y,z), (2.3) 

where Cx is a constant wave speed in the streamwise direction. ChannelFlow seeks solu¬ 
tions of a more general case: 

tr/*^ (m) — It = 0. (2.4) 

Here is the time-ti forward time integration of the Navier-Stokes equations computed 
by a direct numerical simulation (DNS), i.e. f*^{u{t)) = u{t + ti) and cr is a symmetry 
operator to the flow field such that 

a[u, V, w]{x, y, z) = [sxU, SyV, Szw]{sxX + axLx,Syy, s^z + a^^L^). (2.5) 

Here we are following notations for flow symmetries introduced by Gibson et al. (2008). 
The symmetry operator <j consists of two sets of parameters: Sx,Sy,Sz for rotation- 
reflection symmetries (values are either 1 or -1) and ax,az for streamwise and spanwise 
translations (values are real). The symmetry operator a in (2.4) describes the translation 
symmetry of the travelling wave solution after time H. To compute travelling wave so¬ 
lutions propagating in the streamwise direction, the only unknown symmetry parameter 
is the streamwise shift ax{= Cxti/Lx), because, the spanwise shift Oz is set to zero, the 
other symmetry parameters are inherent to the solution and the time shift ti is chosen 
a priori. The parameter Ox is determined as part of the solution process. 

More generally, the symmetries of fluid states can be expressed with the symmetry 
operator (2.5). That is, u = au for certain values of symmetry parameters. The symmetry 
operator a is then expressed in different characters to describe different symmetries of 
fluid states: r for the spatial phase shifts, a for the reflections and s for the shift-reflection 


12 


Jae Sung Park and Michael D. Graham 


or shift-rotation. The four flow symmetries that arise in the present study are: 


ay[u, V, w]{x, y, z) = [u,-v, ■u;](a;, -y, z), 


( 2 . 6 ) 


a:,[u,v,w]{x,y,z) = [u,v,-w]{x,y,-z), 
Txz[u,v,w]{x,y,z) = [u,v,w]{x + ^,y,z + ^) 
si[u,v,w]{x,y,z) = [u,v,-w]{x+ ^,y,-z). 


(2.7) 


( 2 . 8 ) 


(2.9) 


The ay and tr^ symmetries correspond to reflections with respect to the midplanes in the y 
and z directions, respectively. The t^z and si symmetries denote half-domain translations 
in the x and z directions and a shift-reflection symmetry, respectively. In particular, the 
Si symmetry is related to the sinusoidal instability of streaks (Waleffe 1997), which is 
also called the fundamental sinuous mode. 

Finding solutions to equation (2.4) requires good initial guesses. We generate these 
using instantaneous velocity fields from DNS of turbulent trajectories that have been 
symmetrized with respect to the midplane of the domain, y = 0 (i.e. all initial guesses 
satisfy u — Uyu). In particular, since hibernating turbulence has been hypothesized to 
be closely related to travelling wave solutions (Xi & Graham 2012a), we chose initial 
guesses from instants with a lower wall shear stress than the mean value. The time shift 
ti is arbitrary. A relatively large value provides substantial improvement to the rate of 
convergence of the Krylov subspace methods that are used in our computation, but larger 
values of tx require longer to compute . We chose ti = 20 in that it seems to balance 
these two aspects well. The initial guess for the streamwise shift is determined by 
approximating the wave speed as the bulk velocity of the symmetrized initial velocity 
field. Appropriate symmetries are enforced during the time-ti time integration and the 
search procedure. To solve equation (2.4) a Krylov subspace method is used to solve 
the linear systems arising at each Newton step. For better convergence, a trust-region 
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TW 

RCc 

RCt 

Lx 

Lz 

Cx 

Symmetries 

PI 

1315.96 

58.95 

TT 

TT 

0.73 

t *^25 ‘Txz 

P2 

1270.98 

59.88 

TT 

TT 

0.71 

Oy 

P3 

682.57 

43.20 

271 

TT 

0.62 

^ y 5 ‘T'xz ; Si 

P4 

1400 

67.32 

TT 

7r/2 

0.75 

fjy , CTz 

P4SB 

3070 

88.70 

271 

TT 

0.78 

^y t ^ zt '^xz 

P5 

2744.95 

99.02 

TT 

7r/2 

0.78 

(Ty , Cz 


Table 1. Scales and symmetries of travelling wave solutions at their saddle-node bifurcation 
points. The wave speed is normalized by laminar centerline velocity. Note that the bifurcation 
points for P2 and P4 subharmonic branch (SB) correspond to its minimum due to the discovery 
of only one branch. 


limitation to the magnitude of the Newton steps or a hook step within a Krylov subspace 
is computed for the optimal Newton step. The Newton iteration is repeated until an 
accuracy of 0(10“^®) is reached, where the accuracy is the residual of || af*^{u) — w ||, 
using the L 2 norm 


9 11 = 


pi pLx 


2Z/a;Z/2 Jq J _l Jq 


g • g dxdydz 


1/2 


( 2 . 10 ) 


Note that the ChannelFlow code calculates this residual. A detailed description for the 
solution algorithm and DNS can be found in Gibson et al. (2008, 2009) and Viswanath 
(2007, 2009). 
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Figure 2. (Colour online) (a) A bifurcation diagram for five families of travelling wave solutions 
in terms of Reynolds number Rcc = Uch/v and friction Reynolds number Rct = Urhlv. The 
curves for Newtonian turbulence and laminar flows are shown. For P4, solid symbols correspond 
to solutions computed by a subharmonic bifurcation (SB), (b) A Prandtl-von Karman plot for 
a bifurcation diagram. The average velocities as a function of friction Reynolds number are 
shown along with curves for Newtonian turbulence, laminar flow and the Virk maximum drag 
reduction, (c) Only P4 solution family is shown, (d) A comparison of PI and P3 to earlier 
solutions of Waleffe (2001), Nagata & Deguchi (2013) and Gibson & Brand (2014) in the same 


geometry. 
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3. Results and discussion 

3.1. Travelling wave families: bifurcation diagram, mean profiles and overall structures 
We computed five families of nonlinear travelling wave solutions in plane Poiseuille flow, 
which we have labelled PI through P5, in domains of three different sizes: [tt, 2,7r], 
[27r, 2,7r], and [tt, 2,7r/2]. For the same flow geometry, the previous studies of Nagata 
& Deguchi (2013) and Gibson & Brand (2014) used the domain of [ 27 r, 2, tt] and Toh & 
Itano (2003) used the domain of [tt, 2,0.47r]. The minimum spanwise domain size used 
here is about 955,^, corresponding to the length scales of the optimal spanwise wave¬ 
length for ECS (Waleffe 2003) and of the near-wall streak spacing of about 100(5,^ (Smith 
& Metzler 1983). This minimum spanwise length scale is also within the range of the 
critical channel widths about 85-110 wall units used for the minimal flow unit studies 
(Jimenez & Moin 1991). Table 1 presents scales and symmetries of the solutions at their 
bifurcation points. Since only one solution branch is found for P2, the lowest Re solu¬ 
tion is presented. Because we imposed the Oy symmetry on initial guesses, all solutions 
exhibit this ay symmetry; P2 has only this symmetry. The half-period translations in 
the streamwise and spanwise directions, Txz, are found for PI and P3. The shift-reflect 
symmetry si responsible for the fundamental sinuous mode is found for P3. 

The bifurcation diagram for these solution families is shown in figure 2(a). The solutions 
are plotted using friction Reynolds number vs. laminar equivalent Reynolds number. For 
each solution family, a solution with higher i?e,- is an upper branch solution corresponding 
to high drag, while its counterpart is a lower branch solution. For comparison, Newtonian 
turbulence and laminar flow are also drawn. Another representation of the bifurcation 
diagram, a Prandtl-von Karman plot, is shown in figure 2(b). This form is often used 
to represent drag-reduction characteristics in wall-bounded turbulent flows. The bulk 
velocities are plotted as a function of friction Reynolds number along with curves 
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for Newtonian turbulence, laminar flow and Virk MDR. The curve for the Virk MDR 
is generated using its universal mean velocity profile (Virk 1975). We elaborate below 
on the solutions with respect to the Prandtl-von Karman plot. In this representation, 
a “lower branch” solution is above the “upper branch”, because the former has higher 
bulk velocity for the same wall shear stress than the latter. With the exception of the 
lower branch solutions of PI and P3, the maximum Re at which a solution is shown on 
the bifurcation diagram represents the highest Re at which a converged solution could 
be found. Obtaining solutions at higher Re will require further refinements in techniques 
for solving equation (2.4). 

The lower branches of PI and P3 become very close and parallel to the laminar solution 
as Reynolds number increases. Their closeness to the laminar state indicates that they 
are very low drag states. These lower branch solutions have very weak spatial variations 
and Reynolds number dependence and have been successfully continued up to Rer ~ 300 
(corresponding to Rec ~ 40000). Regarding the upper branches of the solution families, 
in the range where we have computed it PI has a similar level of drag (i.e. a similar 
bulk velocity for a given Rer) as Newtonian turbulence. The P3 upper branch, however, 
shows higher drag than Newtonian turbulence, displaying the highest drag level among 
the solutions found in the present study. The P2 solution branch appears to bifurcate off 
PI, a result that is confirmed below when we see that P2 has a broken Uz symmetry. P5 
forms a closed loop (isola). 

Let us now focus on the P4 solution family, which shows very intriguing behaviour 
with regard to Newtonian and viscoelastic turbulence. Figure 2(c) shows only the P4 
solution family. Consider the upper branch (low velocity) solution at the upper range 
of convergence Rcr « 130. This solution branch has a mean velocity very close to that 
of Newtonian turbulence (the black dashed curve). As Re decreases along this branch. 
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the solution remains close to the Newtonian turbulence curve until undergoing a turning 

point at Rcr = 67.32, beyond which the lower branch (LB) appears to approach the Virk 

MDR curve. This solution branch turns around again at Rcr = 88.7, forming another 

lower branch solution (which we call LB2). As we decrease Reynolds number, the new 

branch reaches another bifurcation point at Rcr = 55.63 and we denote the solution 

beyond this point as UB2. Thus, there are three bifurcation points for P4 solution family. 

Interestingly, the bulk velocities at the bifurcation points at Re^- = 55.63 and 67.32 

are remarkably close to the Newtonian turbulence value, while the third bifurcation 

point at Re-r = 88.7 is close to the Virk MDR value. Finally, above the turning point 

at Rct = 88.7, a subharmonic - spatiotemporal period-doubling - branch (SB) arises, 

which has doubled fundamental spatial periods in the x and z directions compared to 

the P4 solution family (i.e. becomes 27r and becomes tt), while the wave speed 

remains constant. Thus at any point in the domain the temporal period of the velocity, 

as measured in the laboratory frame, doubles. The subharmonic solutions are indicated 

by solid symbols. This solution closely follows the MDR curve until Rcr « 105, deviating 

from and then returning to it as Re increases further. 

Prior to proceeding to figure 2(d), it is worth mentioning the linear stability of the solu¬ 
tions. The leading eigenvalues of the solutions are computed in their symmetric subspace 
with Arnoldi iteration (Viswanath 2007). The PI, P3, and P5 lower branch solutions have 
a single, real unstable eigenvalue, while the P4 lower branch (P4-LB) solution has two 
real unstable eigenvalues. The P4-LB2 solution has three real unstable eigenvalues, and 
the P4 subharmonic branch (P4-SB) solution has three real and three complex conjugate 
unstable eigenvalues. Turning from the P4 lower branch to the P4 upper branch (or from 
P4-LB2 to P4-UB2), one real unstable eigenvalue goes complex immediately in a Takens- 
Bogdanov bifurcation (Guckenheimer & Holmes 1983), at which an eigenvalue associated 
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Figure 3. (Colour online) Wave speeds for the travelling waves as a fnnction of Reynolds number 
in (a) outer units and (b) inner units. In (a), the wave speeds are normalized by the laminar 
centerline velocity and the dashed line represents the laminar bulk velocity. 

with a saddle-node bifurcation collides with another eigenvalue. So just beyond their re¬ 
spective turning points, P4-UB and UB2 have one real and one complex conjugate pair 
and two real and one complex conjugate pair of unstable eigenvalues, respectively. The 
PI, P3, and P5 solutions also experience the Takens-Bogdanov bifurcation after crossing 
lower to upper branch. Interestingly, this behavior has also been observed near turning 
points of pipe flow travelling waves (Mellibovsky & Eckhardt 2011; Pringle et al. 2009) 
and thus seems to be rather generic for travelling waves in shear flows. 

Figure 2(d) compares PI and P3 to earlier TW solutions discovered by Waleffe (2001), 
Nagata & Deguchi (2013) and Gibson & Brand (2014) in the same geometry. The curves 
for Waleffe and Nagata & Deguchi solutions are generated from figure 5 of Nagata & 
Deguchi (2013), where the optimal wavelengths were used. Gibson & Brand’s two so¬ 
lutions, named TWl and TW2 in their paper, have the same wavelengths as P3. The 
solution of Nagata & Deguchi, called MS-S in their paper, possesses ay, a^ and si sym¬ 
metries, and Waleffe solution lacks the a^ symmetry compared to MS-S solution. TWl 
has the same symmetries as PI, while the ay symmetry is lost in TW2. The bifurcation 
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Figure 4. (Colour online) Mean velocity profiles for (a) PI, P2 and P3, (b) P5, and (c) P4, in 
comparison to the log-laws for Newtonian turbulence and Virk MDR. (d) Mean profiles of the 
Reynolds shear stress for P4 solution family. A dot-dashed line is a time-average profile for long 
DNS trajectories. 

points of P3, Waleffe and MS-S solutions are very close each other. As shown, the P3 
upper branch appears to be the highest drag state. Note that Waleffe, MS-S and TWl 
solutions were obtained by continuation from plane Couette to Poiseuille conditions, but 
TW2 was obtained by a similar manner to an edge tracking method using a velocity held 
from DNS as an initial guess. To our knowledge, solutions similar to P4 and P5 have not 
previously been found. 

Figures 3(a) and (b) show the wave speed of the solutions as a function of Reynolds 
number in outer and inner units, respectively. In general, the wave speed follows the same 
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trend as the bulk velocity: in a given solution family, a lower branch solution propagates 
faster than an upper branch solution. In outer units in figure 3(a), the laminar bulk 
velocity is plotted in a dashed line. Most of solutions have a larger wave speed than the 
laminar bulk velocity, indicating that they propagate forward when viewed in a reference 
frame moving at the laminar bulk velocity. The wave speeds of the PI and P3 lower 
branch solutions appears to become constant with increasing Reynolds number, while 
those of their upper branch solutions decrease drastically. When plotted in inner units, 
the wave speed behaviour shows almost the same shape as the bulk velocity plot in figure 
2 . 

We now turn our attention to the mean velocity profiles {7+(j/+). Figure 4(a) shows 
these for PI, P2 and P3. For comparison, we also plot the von Karman log-law = 

2.5 lny+-1-5.5 profile of Newtonian turbulence and the Virk log-law = 11.7 ln?/+ — 

17.0 that approximates the mean velocity profile in the MDR regime. As expected from 
the average velocity results in figure 2(b), the lower branch velocity profiles for PI and 
P3 are well above the Virk MDR profile for Rcr > 80 and very close to the parabolic 
laminar profile. The highest-drag solution of the P3 upper branch solution shows a mean 
velocity profile well below the von Karman log-law profile. The velocity profile for P2 
shows a similar character to the PI and P3 upper branch profiles. 

In figure 4(b), the mean velocity profiles for P5 are shown at its minimum and max¬ 
imum Reynolds numbers, and at Rcr « 104. In particular, the lower branch velocity 
profile at Rcr = 104.51 very closely approaches the Virk MDR log-law in the range 
15 < 2/+ < 45. 

Figure 4(c) shows mean velocity profiles for P4. Starting from its first bifurcation 
point at Rcr = 67.32, the upper and lower branch profiles seem to approach toward 
two distinct limits, the von Karman and Virk MDR profiles, respectively, as Reynolds 
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number is increased. In particular, the lower branch profile at Rcr = 88.38 very closely 

approaches the Virk MDR log-law profile over a relatively wide range, 15 < y’*' < 70. 

This lower branch profile is very similar to the conditionally-sampled DNS velocity profile 

and experimentally observed profile for Newtonian hibernating turbulence (Xi & Graham 

20126; Whalley et al. 2014). Meanwhile, the upper branch profile lies very close to the 

experimentally-observed mean profile of Newtonian turbulent flows and approaches the 

von Karman log law profile at large as does active turbulence. The subharmonic 

branch profile is also presented for Rcr = 100.66, showing a similar shape as the lower 

branch. Hence, the upper, lower and subharmonic branches of the P4 solution family 

may represent an envelope in state space encompassed by the mean profiles of both the 

Virk MDR and the classical Newtonian turbulence. 

For the P4 solution family, we plot Reynolds shear stress profiles in figure 4(d), in 
comparison to the time-average Newtonian profile for long DNS trajectories. As Reynold 
number is increased, the Reynolds shear stress for lower branch solutions decreases, while 
the upper branch profile increases. Compared with the time-average Newtonian profile 
at Rcr ~ 86, the reduction in the lower branch solution is substantial. Interestingly, 
the upper branch profile is slightly higher compared with magnitudes of the Newtonian 
profile in y+ < 30, but it almost collapses onto the Newtonian profile for y+ > 45. 

Now we examine the streak structures, as represented by the contours of streamwise 
velocity fluctuations in the x — z plane at y = —0.5; these are shown in figure 5. Except 
for P2, lower branch solutions are presented. The low-speed and high-speed streaks are 
denoted by negative (blue) and positive (red) fluctuations, respectively. A subharmonic 
sinucose mode (Waleffe 1997), which is a combination of sinusoidal for the low-speed 
streak and a varicose mode for the high-speed streak, is identified for PI, P3 and P5. A 
fundamental sinuous mode is observed for the P2 solution (figure 5(b)), while P4 exhibits 



Figure 5. (Colour online) Contours of streamwise velocity fluctuation in the x — z plane at 


y = —0.5 for (a) PI, (b) P2, (c) P3, (d) P4 at Rec = 1800, and (e) P5 at Rec = 3600. Lower 


branch solutions are presented except for P2, for which there is only one branch. Blue and red 


indicate low- and high-speed streaks, respectively. 


a fundamental varicose mode (figure 5(d)). From a flow symmetry point of view, the 


symmetry is clearly seen for PI, P4 and P5, whereas this symmetry is broken for P2 and 


P3. The Txz symmetry for PI and P3 is also identified in hgures 5(a) and (c). 


To clearly illustrate the subharmonic bifurcation arising on the P4 lower branch around 
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Figure 6, (Colour online) (a) P4 lower branch solution at Re^ = 2750 (two periods in x and 


a are shown) and (b) subharmonic solution at Rcc = 3800, Shown are colour contours of the 


wall-normal velocity in the x — z plane at y = —0,5, The solid-outlined box at lower left shows 


the size of the domain in which the P4 solutions are found. The subharmonic solution has a 


broken mirror-symmetry about the midplane of this box in the 2 direction (dashed line), as 


well as a broken discrete translation symmetry in both the x and 2 directions - the solution in 


TT < a: < 2?! is shifted by tt in the 2 -direction from the solution in 0 < a: < tt. 


Rcc = 3070, we plot in figures 6 (a) and (b) the wall-normal velocities on the x-z plane 
of solutions at Rcc = 2750 and 3800, respectively. The fundamental spatial periods of 
the P4 lower branch solution are = tt and = 7r/2; figure 6(a) shows two periods 
of this solution in each direction - the unit cell is the solid-outlined box at lower left. In 


the subharmonic solution, figure 6(b), mirror-symmetry with respect to the z direction 
midplane (dashed line) of the unit cell of P4 is broken, even though the symmetry still 
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Figure 7. (Colour online) Vortical structures of the travelling wave solutions as illustrated by 
the swirling strength Xd- (a) lower and (b) upper branch solutions for PI, (c) P2 solution, 
(d) lower and (e) upper branch solutions for P3, (f) lower and (g) upper branch solutions for 
P4, (h) lower and (i) upper branch solutions for P5. P5 is shown at Rcc = 3600, all others at 
Rcc = 1800. The bottom-half channel is only shown due to the mirror-symmetry with respect 
to the channel center. The numbers in parenthesis are the maximum swirling strength. Dark red 
tubes are isosurfaces at 2/3 of the maximum swirling strength, and transparent blue isosurfaces 
indicate critical layer surfaces, where the local streamwise velocity matches the wave speed. 

holds for its own (larger) fundamental domain. In addition, the subharmonic solution has 
a broken discrete (half domain shift) translation symmetry in both x and z directions: 
the solid-outlined box is not same as the solution in tt < a; < 27r and 0<z<7r/2orin 
0 < a; < TT and 7r/2 < z < tt. 

The streak structure of a flow is closely related to the streamwise vortical structure. 


Figure 7 shows contours of swirling strength Ad, the imaginary part of the complex 
conjugate eigenvalues of the velocity gradient tensor (Zhou et al. 1999), in the bottom 
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half of the channel. The contours represent 2/3 of the maximum swirling strength for 

each solution, which is presented in parentheses in the figure. In a given solution family, 

lower branch solutions have weaker vortex strength than upper branch solutions. We also 

depict the critical layer surface, where the local steamwise velocity equals the wave speed, 

u{x,y,z,t) = Cx (Maslowe 1986; Hall & Sherwin 2010). The PI and P3 lower branch 

solutions, which have the subharmonic sinucose mode, show similar vortical structures. 

The vortex cores appear to expand between just above critical layer and the channel 

center, whereas the vortex cores of the PI and P3 upper branch solutions are located 

very close to the critical layer. The fundamental sinuous mode of P2 displays staggered 

vortices, which are also located close to the critical layer. The fundamental varicose 

mode of the P4 lower and upper branch solutions displays different inclination angles of 

the vortex legs with respect to the wall: the vortex legs of the lower and upper branch 

solutions are inclined about 20 and 40 degrees to the wall, respectively. The vortex 

cores are also located around the critical layer. This vortical structure - which has the 

same symmetry as a hairpin but does not display a “head” - is also observed for other 

travelling wave solutions in the same geometry (Gibson & Brand 2014) and plane Couette 

flow (Itano & Generalis 2009; Deguchi & Nagata 2010). 

3.2. Travelling wave structure and critical layers 

As described in the Introduction, prior work has addressed the structure and mechanism 
of nonlinear travelling waves in the context of nonlinear critical layer dynamics and in 
particular has noted the role that streamwise-wavy structures localized near the critical 
layer play in the self-sustaining process of at least one family of ECS (Wang et al. 2007; 
Hall & Sherwin 2010). Therefore it is of interest to illustrate the channel flow ECS found 
here in relation to the critical layer position. In order to do so, in figure 8 we calculate 
velocity deviations u{x = Xc,y,z,0) — U{y,z) in the y — z plane of the PI, P3 and 
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Figure 8. (Colour online) Contours of streamwise velocity deviations u(x = Xc , y, z,t) — U{y, z) 
in the y — z plane, where U(y, z) is the streamwise-averaged streamwise velocity. The streamwise 
location x^ is chosen to illustrate the distinguishing featnres of the deviation. The flow fields 
shown are (a) lower and (b) upper branch solutions for PI, (c) lower and (d) upper branch 
solutions for P3, (e) lower, (f) upper and (g) subharmonic branch solutions for P4. Except for 
the subharmonic solution, which is shown at Rcc = 3600, all solutions are at Re^ = 1800. The 
black line represents the critical layer in the y — z plane a.t x = Xc- 

P4 lower and upper branch solutions along with the critical layer position (black thick 
line), where u{x = Xc,y,z,0) = Cx- Here U{y,z) is the streamwise-averaged streamwise 
velocity {not the streamwise- and spanwise-averaged velocity Um{y))- As the velocity 
deviation varies along the streamwise direction, a location Xc in the streamwise direction 
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is chosen so that the distinguishing features of the deviation are best illustrated. The 

full time-dependence of these structures is shown in supplementary movies 1-7. The 

PI and P3 lower branch solutions exhibit relatively strong deviations near the channel 

center as shown in figures 8(a) and (c), (and supplementary movies 1 and 3) from which 

observation we may call them “core modes”. In particular, the P3 solution shows well- 

localized deviations near the channel center. Upper branch solutions exhibit stronger 

deviations throughout the channel height compared to lower branch solutions. In both 

cases, however, the fluctuations seem to be bounded between the top and bottom critical 

layers (also see accompanying online movies 2 and 4 in the supplementary materials for 

figures 8(b) Pl-UB and (d) P3-UB, respectively). 

P4 shows a different structure. Figures 8(e), (f) and (g) (and supplemental movies 
5-7) show the streamwise velocity deviations for its lower and upper branch solutions at 
Rbc = 1800 and for a subharmonic solution at Rcc = 3600, respectively. The deviations of 
the lower and subharmonic branch solutions are highly localized very close to the critical 
layer, consistent with the Couette flow ECS results of Wang et al. (2007) and Hall & 
Sherwin (2010). For the upper branch solutions, while strong deviations are observed 
across much of the channel, they are clearly organized by the critical layer. The clear 
organization of deviations around the critical layer suggests a connection to the critical 
layer dynamics, based on which this solution family may be called a “critical layer mode”. 
As seen in the vortical structures in figures 7(h) and (i), the structure of P5 is also 
strongest around the critical layer. 

This distinction between core modes and critical layer modes does not seem to have 
been previously identified. The lower branch Couette (Blackburn et al. 2013; Hall & 
Sherwin 2010; Wang et al. 2007) and pipe (Viswanath 2009) flow travelling waves stud¬ 
ied previously would be classified as “critical layer” rather than “core” modes, and the 
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mechanistic studies of Blackburn et al. (2013) and Hall & Sherwin (2010) are focused on 
critical layer modes. 

In addition to movies for streamwise velocity deviations, there are additional accom¬ 
panying online movies 8-17 in the supplementary material for the PI - P5 travelling wave 
structures. These are shown in the y — z plane, where streamwise velocity is represented 
by colour contours, wall-normal and spanwise velocities are shown by arrows and the 
critical layer is shown as a black curve. 

3.3. Exact coherent states on the laminar-turhulent boundary 
Some of the lower branch travelling wave solutions in wall-bounded turbulent flows are 
embedded in the laminar-turbulent boundary (Skufca et al. 2006; Wang et al. 2007). If 
such solutions have only one unstable direction (Skufca et al. 2006), which indicates that 
they are stable with respect to perturbations on the boundary, they are edge states. By 
combining linear stability analysis of the nonlinear travelling wave families with direct 
time-integration of initial conditions perturbed along unstable directions of the travelling 
waves, we have determined whether the solutions live on the basic boundary and whether 
they are edge states. Specifically, an initial condition for a trajectory is generated by ad¬ 
dition of a small perturbation along an unstable eigenfunction to a lower branch solution. 
Both positive and negative perturbations are considered, and if there is an unstable di¬ 
rection for which the trajectory starting on one side of the lower branch solution develops 
into turbulence, while the other decays directly to the laminar state, then the solution 
lives on the basin boundary. If, additionally, there is only one unstable eigenvalue, then 
the travelling wave is attracting in all other directions besides its unstable one and is 
thus an edge state. 

Figure 9(a) shows dynamical trajectories along the unstable directions for the P3 and 
P4 lower branch solutions at Rec = 1800 projected onto the plane of energy dissipation 
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Figure 9. (Colour online) (a) Time evolutions in the energy input rate and dissipation rate for 
DNS trajectories starting from nonlinear travelling waves P3 (o) and P4 (o) perturbed along an 
unstable eigendirection at Rcc = 1800. The solid and open symbols correspond to lower branch 
and upper branch states, respectively. Along trajectories, dot spacing is At = 2. The laminar 
state is at (1,1). The dashed line represents D = I. (b) Bifurcation diagram with solutions on 
the basin boundary shown in colour. 

rate (D) and energy input rate (7), 

J {pu\x=o - pu\x=Ljdydz. (3.2) 

Both values are normalized by their laminar values such that the laminar state is at 
(1,1). Recall that the fundamental domain size is different between P3 and P4: = 2tt 

and Lz = TT for P3 and Lx = tt and Lz = 7r/2 for P4. In both cases, one perturbation (the 
one that increases I and D) leads to turbulence and the other to laminar, indicating that 
these travelling waves are on the basin boundary for turbulent flow in their respective 
domains. Even though there are multiple unstable eigenvalues for the P4 lower and 
subharmonic branches, only one (real) unstable eigenvalue gives the aforementioned two 
types of trajectories of perturbations. Figure 9(b) shows the results of this analysis for the 
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entire bifurcation diagram - the solutions shown in gray are not on the basin boundary 
while other others are - we find that only P3 and P4 display parameter ranges where 
they are on the boundary. For P3 this range is 50 < Rcr ^ 123 (1200 < Rcc ^ 5000). For 
P4, both lower and subharmonic branch solutions are found to lie on the basin boundary 
for their respective domains. The ranges are 69 < Rcr < 88.7 (1600 < Rcc < 3070) 
for P4-LB, 61 < Rcr < 88.7 (1400 < Rcc < 3070) for P4-LB2, and 88.9 < Rcr < 
107 (3080 < RCc < 4000) for P4-UB. Furthermore, since it is found that the P3 lower 
branch solutions have only one unstable eigenvalue in the symmetric subspace, they are 
indeed edge states in that space. Even though there is only one eigenfunction that gives 
rise to an escape scenario from the basin boundary for the P4 lower branch solutions, 
they have multiple unstable eigenvalues so are not edge states. 

3.4. Connections between travelling wave solutions and turbulent trajectories 
One motivation for studying nonlinear travelling waves in shear flows is the idea that 
these states form the state-space skeleton of the turbulent dynamics. In this section, we 
address this issue, examining how close the turbulent trajectories approach the travelling 
wave solutions. We focus on the domain x Ly x Lz = x 2 x 7r/2, the same box size 
as the P4 solution family. We choose Rcc = 1800 {Rer = 85) and perform simulations 
at constant mass flux. In contrast to the travelling wave computations, these turbulence 
simulations are performed without imposing any symmetries on the flow. Comparisons 
are made with travelling waves that have the same Re^- Two methods of comparison are 
used: the hrst examines the probability distribution function for the instantaneous mean 
velocity profile and the second projects the state space dynamics into three dimensions 
corresponding to physically meaningful averaged quantities. 

Figures 10(a) and (b) show the probability density functions (PDFs), plotted on a 
logarithmic scale, for the mean velocity prohles at each wall normal position y (figure 
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Figure 10. (Colour online) The probability density function of mean velocity profile from 
DNS, along with P4 lower branch (P4-LB), another lower branch (P4-LB2), and upper branch 
(P4-UB) solutions on (a) outer scaling and (b) scaling (inner nnits based on instantaneous 
area-averaged wall shear stress). The black line is the time-averaged DNS velocity profile. A 
logarithmic scale is used, with blue indicating vanishing probability. 

10(a)) or y* (figure 10(b)) in DNS based on outer units and ‘^’-scaling (instantaneous 
inner scales), respectively. The PDF is normalized so that the integral over the whole 
PDF equals 1. Using outer scaling (non-time-dependent scaling), it is difficult to com¬ 
pare a DNS trajectory to TW solutions because each TW has different friction velocity. 
However, as highlighted by previous studies (Xi & Graham 20126; Agostini & Leschziner 
2014), the ‘^’-scaling, which leads all profiles to collapse to the same curve near the wall, 
is the proper one to use for instantaneous quantities with which a TW solution can be 
directly compared to an instantaneous flow field. Here, we used DNS results for 40000 
time units to compute PDFs. According to Xi & Graham (20126) and our calculations, 
approximately 8 - 9% of the total simulation time is spent near the Virk-like state. Thus 
the data is sufficient to capture approaches to the Virk log-law in PDFs. In figure 10(a) we 
see that near the wall, the DNS velocity profile is very nearly bracketed between the P4 
upper branch and lower branch solutions, while deviations from these solutions become 







32 Jae Sung Park and Michael D. Graham 

more prevalent near the center. The same trend is apparent in the plot in instantaneous 


inner units, figure 10(b), which emphasizes the strong similarities in the near-wall be¬ 
haviour as well as the transient approaches of the DNS mean velocity profile toward the 
P4 lower branch solutions. It appears that for y* < 30 the P4 travelling waves form 
an approximate envelope for the PDF of the DNS mean velocity profile. Furthermore, 
relatively high probability regions (red) are observed around the von Karman log-law 
and P4-UB solution. Interestingly, there is also a slightly high probability region (yellow) 
close to the P4-LB solutions. 

Now we visualize the approach of turbulent trajectories to travelling wave solutions in 
state space. To do so, we project turbulent trajectories onto a three-dimensional space 
using the following quantities: disturbance kinetic energy (KE), energy dissipation rate 
(D), and area-averaged instantaneous wall shear stress normalized by its mean value 
(tw/tw)- The disturbance kinetic energy is defined as follows: 



(3.3) 


where uiam is the parabolic laminar profile. Figure 11 shows a turbulent trajectory as 
well as the P4 travelling waves projected onto these three quantities, as well as the 
joint probability density function (PDF) of KE and D at the bottom of the figure. 
Note that all quantities are calculated for only the bottom half of the domain. The 
dynamical trajectory spends most of its time within one core region of state space, which 
we can identify with normal or “active” turbulence; the P4-UB solution is in this region, 
as also seen on the joint PDE. The trajectory occasionally escapes, however, from the 
active region, approaching the P4 lower branch solutions. During these excursions, some 
trajectories pass through the vicinity of P4-LB, approaching P4-LB2 very closely. When 
returning to the active region, the path strongly overshoots the core of these region 


and might be considered as a turbulent burst. Similar observations on the relationship 



Exact coherent states in Poiseuille flow 


33 


1 . 2 > 



KE 


Figure 11. (Colour online) A state-space visualization of DNS trajectories, projected onto three 
dimensions: disturbance kinetic energy (KE), an energy dissipation rate (D), and normalized 
instantaneous wall shear stress (tto/t™). The grey line indicates the turbulent trajectory, to 
which black dots are attached at intervals of Ih/Uc- A joint probability function of KE and D is 
shown at the bottom of the figure. The labelled symbols (<) are P4 solutions. Points (i), (ii), and 
(iii) are the closest visits to P4-LB, P4-LB2, and P4-UB, respectively. Note that all quantities 
are calculated only for the bottom half of the channel. 

between invariant states and bursting events have been made for Couette flow (Kawahara 
& Kida 2001) and channel flow (Toh & Itano 2003). The closest visits to P4-LB, P4- 
LB2, and P4-UB, respectively, are labelled as points (i), (ii), and (iii). The mean velocity 
profiles at these three instants and for the P4 travelling waves are plotted in figure 
12(a). The profiles for instants (i) and (ii) appear very similar to the P4-LB and P4-LB2 
solutions, respectively, while instant (iii) has a similar proHle to the P4 upper branch 
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Figure 12. (Colour online) (a) Instantaneous mean velocity profiles for instants (i), (ii), and 
(iii), where the closest approach to P4 lower branch (LB), second lower branch (LB2), and upper 
branch (UB) solutions are observed, respectively, (b)-(c) Flow structures at instants (i) and (iii). 


and the von Karman mean profile. Flow structures for instants (i) and (iii) are visualized 
in figures 12(b) and (c), where to facilitate comparison we use the same vortex strength 
and critical layer isosurfaces as in figures 7(f) P4-LB and (g) P4-UB. Note that flow 
structure for instant (ii) is very similar with that of instant (i), but shows weaker vortex 
motion. There is substantial similarity between the vortical and critical layer structures 
of the snapshots and the travelling waves. The critical layer for instant (i) shows weak 
streamwise variation and the vortex motions are seen to be localized around the critical 
layer as also seen for P4-LB. Similarly, instant (iii) and P4-UB resemble one another, 
particularly with regard to their inclined vortical structures. 

To further address the relationship between the DNS trajectories and travelling waves, 
we calculated a norm of the difference between a DNS velocity field u and a travelling 
wave. Taking into account the arbitrary phase in x and z of the velocity fields, this 
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Figure 13. (Colour online) Distances <5 between a DNS trajectory and P4 travelling wave solu¬ 
tions. Instants (i), (ii), and (iii) are the closest visits to P4-LB, P4-LB2, and P4-UB, respectively, 
indicated on figure 11. 

distance, which we denote 6, is calculated as follows: 


6{t) = min min \\ u{x + xo,y, z + zo,t) - UTw{x,y, z) \\, (3.4) 

where xq and zq are phase shifts in the x and z directions, respectively. Time series of S{t) 
computed for the P4 solution family are shown in figure 13. The distances between P4- 
LB and instant (i), P4-LB2 and instant (ii), and P4-UB and instant (iii) are 5.7 x 10“^, 
6.2 X 10“^, and 3.5 x 10“^, respectively. Those values are comparable to distances at 
the points of the closest visits to travelling waves for pipe flow (Viswanath & Cvitanovic 
2009) and to equilibria for Couette flow (Halcrow et al. 2009), which are of order 0(10“^). 
Thus, the closeness of turbulent trajectories to P4 travelling waves is identified using full 
velocity fields for minimal channel flow. 

The state-space picture that clearly shows close approaches to multiple TW solutions 
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and the closeness to the TW solutions using full velocity fields have yet to be reported 
in the channel flow literature. Furthermore, it must be emphasized that those multiple 
travelling waves belong to the same solution family. 

These results confirm the hypothesis posed in prior work (Graham 2014; Xi & Gra¬ 
ham 20106 , 0 , 20126) that the “active” and “hibernating” phases of minimal channel 
turbulence correspond to time intervals where the trajectory is close to upper and lower 
branch travelling waves, respectively. Finally, returning to the bifurcation diagram we 
recall that there are also two upper branch P4 solutions, but we were only successful in 
computing one of them over a broad range of Rcc- We might speculate that if the second 
upper branch solution could be found at the Rcc at which we have performed the DNS 
it would also lie in the core active turbulence region. 

4. Conclusion 

We have computed five new families of nonlinear travelling wave solutions, denoted 
P1-P5, in Newtonian plane Poiseuille flow. As Reynolds number is increased, the PI 
and P3 lower branch solutions become close and parallel to the laminar solution branch, 
indicating that they are very low-drag states. The P2 solution branch results from a 
symmetry-breaking bifurcation from PI. P5 forms a closed loop (isola). Most interest¬ 
ingly, the P4 solution family shows very intriguing behaviour in terms of mean properties 
as Reynolds number is increased. The average velocities of the lower and upper branches 
appear to approach the Virk MDR profile observed in viscoelastic turbulence and the 
classical Newtonian (von Karman) profiles, respectively. The former observation adds to 
the set of results in which mean velocity profiles close to the Virk profile are found in 
Newtonian flow (Kerswell et al. 2003; Bandyopadhyay 2006; Xi & Graham 20106, a; Du- 
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bief et al. 2011; Xi & Graham 2012a, 6). On the lower branch, a subharmonic bifurcation 

arises around Rcr ~ 90, giving rise to spatiotemporal period doubling. 

The structures and symmetries of the various solution families are described. The fluc¬ 
tuations of the PI and P3 solutions are largest near the channel center so we have denoted 
them as core modes, while the P4 and P5 solutions display fluctuations localized around 
the critical layer so we call them critical layer modes. Over a range of Reynolds numbers 
P3 and P4 lower branch solutions are embedded in the laminar-turbulent boundary. 

Finally, we addressed the issue of how close the turbulent trajectories approach the 
travelling wave solutions, focusing on the P4 family. In prior work (Graham 2014; Xi & 
Graham 20106, a, 20126) it was hypothesized that “active” and “hibernating” phases of 
minimal channel turbulence correspond to time intervals where the trajectory is close to 
upper and lower branch travelling waves, respectively. The present results corroborate 
this hypothesis. The turbulent trajectory spends most of its time within a region of state 
space that can be identified as normal or “active” turbulence; the P4-UB solution is 
in this region, while the hibernating intervals are approaches to the P4 lower branch 
solutions. 
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